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We investigate theoretically the stationary states of two bosons in a one-dimensional optical lattice 
within the Bose-Hubbard model. Starting from a finite lattice with periodic boundary conditions, we 
effect a partial separation of the center-of-mass and relative motions of the two-atom lattice dimer 
in the lattice momentum representation, and carefully analyze the eigenstates of the relative motion. 
In the limit when the lattice becomes infinitely long, we find closed-form analytic expressions for 
both the bound state and the dissociated states of the lattice dimer. We outline the corresponding 
analysis in the position representation. The results are used to discuss three ways to detect the 
dimer: by measuring the momentum distribution of the atoms, by finding the size of the molecule 
with measurements of atom number correlations at two lattice sites, and by dissociating a bound 
state of the lattice dimer with a modulation of the lattice depth. 
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I. INTRODUCTION 

Optical lattices containing quantum degenerate Bose 
and Fermi gases [1] have been a major topic in atomic, 
molecular, and optical physics of late, one motivation 
being experimental realizations of long-standing lattice 
models in condensed-matter physics and statistical me- 
chanics [2, 3]. In an optical lattice the transition ampli- 
tude (hence, probability) for an atom to tunnel from one 
lattice site to the next can be tuned over a wide range 
by adjusting the intensity of the lattice light. Feshbach 
resonances also permit a precise broad-range adjustment 
of the atom-atom interactions [4, 5]. This means that 
both the effective mass of the atoms and the atom-atom 
interactions are subject to experimental control. 

Our broad theme is aggregates of atoms-call them lat- 
tice molecules-that can be formed and controlled in an 
optical lattice by controlling both the one-particle prop- 
erties and the atom-atom interactions. The simplest one 
is the composite of two atoms, which we and others term 
a "lattice dimer." Literature on variations of the lat- 
tice dimer is accumulating steadily (see [6-17] and refer- 
ences therein) and the trimer problem has also been ad- 
dressed [18]. It seems to us, however, that in the general 
quest toward novel physics seemingly elementary aspects 
such as the stationary states of the dimer and their limits 
when the lattice becomes infinitely long have been over- 
looked. The purpose of the present article is to remedy 
the situation. 

We analyze the eigenstates of the Bose-Hubbard 
Hamiltonian for two bosonic particles in a lattice. We 
start in Sec. II with a finite lattice, employing periodic 
boundary conditions. We transform the Hamiltonian into 
lattice momentum representation and separate what are 
the lattice analogs of the center-of-mass motion and rel- 
ative motion of the atoms in a two-atom lattice dimer. 
The separation is, however, not complete [7, 9, 14]: In 
the analysis of the relative motion, the amplitude of the 
tunneling between lattice sites will get scaled by a fac- 
tor that depends on the center-of-mass motion. We, of 
course, have a numerical solution of the time-independent 



Schrodinger equation on hand, but also discuss qualita- 
tive features of the solution analytically. One bound state 
is found regardless of whether the atom-atom interac- 
tions are attractive or repulsive [8], along with what is 
the finite-lattice analog of the dissociation continuum of 
the lattice dimer. A surprising number of mathematical 
issues arise from the discreteness of the lattice and the 
boson symmetry of the states. We clarify them carefully. 

The limit of a long lattice is the subject of Sec. HI. 
The main item here is the development of a closed-form 
analytic expression for the dissociated state of the dimer 
in the lattice momentum representation. Again, in con- 
trast to past discussions in terms of Green's functions 
or scattering theory, we are after stationary states of the 
time-independent Schrodinger equation. Many past pub- 
lications analyze the limit of a long lattice, occasionally 
without ever saying so. Inasmuch as a direct compari- 
son is possible, our results agree with the previous re- 
sults [8, 14]. 

Section IV presents the analysis of both the finite lat- 
tice and the long-lattice limit in the position (lattice site) 
representation [14]. We reaffirm the results found from 
the momentum representation, and also find additional 
results such as an analytical approximation to the en- 
ergies of the would-be dissociated states when there are 
many lattice sites. 

As an example of the utility of our analytical results, 
we discuss possible ways of detecting the lattice dimer 
in Sec. V. Momentum distribution of the atoms [8], pair 
correlations of atom positions, and dissociation rate of 
a bound dimer when the depth of the lattice is modu- 
lated [8] are the specific cases. The remarks in Sec. VI 
conclude the article. 



II. FINITE LATTICE 

Let us start with the Bose-Hubbard model, ostensi- 
bly in one dimension, though the same mathematics ap- 
ply with straightforward modifications in multiple dimen- 
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sions. The Hamiltonian with standard conventions reads 
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The index k runs over the lattice sites, L of them; k = 
0, . . . ,i — 1. We use periodic boundary conditions, so 
that k ~ L is the same as fc = 0, and hkewise fc = — 1 
means k = L — 1. This would be physically valid for a 
ring lattice. We do not intend to restrict the discussion 
to ring lattices, but the motivation is twofold. First, in a 
long-enough lattice the results will obviously be similar 
for arbitrary boundary conditions. Second, we have a 
sound framework to go mathematically to the limit of an 
infinitely long lattice. 

The periodic boundary conditions bring in some subtle 
topology that we will not address but rather hide with 
our subsequent choices of the parameters. The principal 
one is that the number of lattice sites L is taken to be 
even. This, again, is something that cannot materially 
influence the results for a very long lattice. The sign of 
the site-to-site hopping amplitude J is a matter of choice 
and can be flipped with the trivial canonical transforma- 
tion Qk — > (— l)'^afe. Here we take the native hopping 
amplitude to be positive, J > 0. A similar freedom does 
not apply to the strength of atom-atom interactions at 
each site; U > and U < correspond to repulsive and 
attractive interactions between the atoms, respectively. 

In order to take advantage of the translation invari- 
ance, we next convert to momentum representation. The 
operators 
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are also boson operators when q runs over a suitable 
set of values, for instance, q = 2ttQ/L with the inte- 
gers Q picked so that we have L quasimomenta in the 
first Brillouin zone of the lattice. For an even num- 
ber of sites L our standard choice is to allow the values 
Q = -L/2 + 1, -L/2 + 2,..., L/2; L of them. The ad- 
dition of quasimomenta q is understood modulo 27r, so 
that the result belongs to the same set. The standard 
discrete Fourier transform relations 

^e'''^-=W[,,o], 5]e^«'==W[,,o], (3) 

k q 

where the brackets remind us of the fact that comparison 
of quasimomenta is modulo 2tt and comparison of lattice 
sites modulo L, give the inverse of the definition (2) 
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The Hamiltonian in momentum representation then be- 
comes 



91,92,93,94 



where the transformation has diagonalized the term in 
the Hamiltonian describing site-to-site hopping, 

ujq = — J cosq . (6) 

We study the most general state vector for two bosons, 



E 

Pl,P2 



A(pi,p2)4^ctjO) 



(7) 



Here pi and p2 again are lattice momenta from the first 
Brillouin zone, and |0) is the particle vacuum. As the 
two creation operators commute, we assume without any 
loss of generality the symmetry of the A coefficients 
A{pi,p2) = A{p2,pi). Action of the atom-atom inter- 
action part in the Hamiltonian then leads to 
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E ^[91- 
91,92,93,94 

= 2 E ^(pi'P2)'5[gi+g2,pi+P2] 4i42 |o) ; (8) 

9l,92,Pl,P2 

see the Appendix. The atom-atom interaction therefore 
conserves the total lattice momentum of the atom pair, 
modulo 27r. 

Suppose the state vector is confined to the subspace 
with the total momentum P and write pi^2 = ± 9- 
For clarification, we note that we refer to pi, p2, and 
q as quasimomenta and P as the total momentum even 
though all of the quantities are obviously dimensionless. 
This notation introduces a set of lattice-momentum-like 
quantities q to describe the relative motion of the two 
atoms in such a way that the state (7) reads 



E^(9)4 



where A{q) now has the symmetry 
A{q) - A{-q) . 



(9) 



(10) 



(5) 



As always, the sums ip ± q are modulo 27r so that they 
have the proper values of a lattice momentum in the first 
Brillouin zone. 

There are two details to consider. First, as Eq. (8) 
shows, the comparison of the total momenta P in atom- 
atom interactions is modulo 27r. Two particles at the 
edge of the Brillouin zone can scatter into two particles 
at the center of the Brillouin zone. On the other hand, 
ip stands for the average momentum of the two atoms, 
and two atoms at the center of the Brillouin zone {P = 0) 
clearly have a different physical signature than two atoms 
at the edge of the Brillouin zone {P = 2tt). The proper 
range of P is therefore {—2tt,2tt]. Second, the newly 
introduced summation index q may have some peculiar 
properties. If ^P in itself is a legal lattice momentum, 
then so is q, and the sum in (9) is taken to run over the 
first Brillouin zone as usual. In contrast, if ^P is not a 
lattice momentum, the summation indices q are of the 
form q ~ 2tiQ / L with a half-integer Q, and the sum is 
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taken to run over the values Q = — L/2 + 1/2, . . . , L/2 — 
1/2. We then say that the relative motion is governed by 
half-integer lattice momenta. 

At least in principle, a state with a half-integer lattice 
momentum may be prepared experimentally. Suppose 
that the atom-atom interactions arc first turned off, for 
instance, by means of a Feshbach resonance, and two 
atoms are put in the states with = 0, q2 = 2tt/L. 
This is a state with a half-integer relative momentum. 
Moreover, if subsequently the atom-atom interactions are 
turned on so that the lattice translation invariance is not 
broken in the process, a nontrivial half-integer state is 
liable to arise for the relative motion of the two atoms. 
In the limit of an infinitely long lattice the difference be- 
tween ordinary and half-integer lattice momenta must be 
irrelevant, but for completeness we occasionally mention 
it. 

Let us continue with an ansatz of the form (9). In 
view of Eq. (8), the action of the Hamiltonian on the 
state vector will give 



two atoms. Using (6) we have 



jP+q ^P-q 



X c 



fP+g ^P-q 



|0) 



(11) 



where q' is a legal or a half-integer lattice momentum 
exactly like q. The time-independent Schrodinger equa- 
tion H\iIj) = E\il)) is obviously satisfied if the coefficients 
of the products of the boson creation operators are the 
same for each q on both sides of the equation, which is 
the case if the coefficients A{q) satisfy 



- 1 P,, + - 1 p_ J + ^ E ^('^O - f ^('^) • (12) 

q' 



Under the restriction (10), Eq. (12) is not only 
sufficient but also necessary for the solution of the 
Schrodinger equation. This could be shown by restricting 
the summation index q in the ansatz (9) to non-negative 
values only. The ambiguity in the argument that the 
states c\ c\ |0) are the same for q and —q would 

thereby be removed, but at the cost of severe complica- 
tions elsewhere. We therefore retain the double counting 
of the states in the formal analysis. So far we have re- 
stricted the values of the summation indices q and q' to 
the first Brillouin zone. However, if convenient, A[q) may 
be regarded as a periodic function of q with the period 
of 2tt. In that case the sums run over appropriate dis- 
crete values within a period, including an end point of 
the period only once. 

Unlike for free atoms with a quadratic dispersion rela- 
tion when the total momentum simply provides an addi- 
tive constant to the energy, here the total momentum- 
center-of-mass momentum if you will-has a profound ef- 
fect on the stationary states of the relative motion of the 
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= — 2Jcos(ip) cosg. 
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so that the center-of-mass motion couples to the relative 
motion described by the quasimomenta q and effectively 
scales the hopping amplitude J. The center-of-mass mo- 
tion and the relative motion do not completely sepa- 
rate [7]. For future use we define the frequency scale of 
the relative motion of the atoms for a given total center- 
of-mass quasimomentum P as 



Op = 2Jcos(iP) . 



(14) 



We also scale the dimensional parameters in the problem 
to fip, defining the dimensionlcss quantities representing 
energy of a state and strength of the atom-atom interac- 
tions as follows: 



Mlp ' f2p 



(15) 



With these ingredients the solution of the time- 
independent Schrodinger equation is reduced to two 
steps. First, Eq. (12) gives 
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and this further gives the equation for the energy eigen- 
value 



I .1 

L ^ E/h + 2Jcoai^P)cosq 
In the dimensionlcss variables (15), Eq. (17) reads 



L ^-^ UJ + cos q tC 

q 
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(18) 



Whether we deal with integer or half-integer quasi- 
momenta, we may always pair them up in the sum in 
Eq. (17) so that for every q there is a q' such that 
cos q ~ — cos q' . This implies that the sign of the fac- 
tor 2Jcos(iP) has no effect on the possible solutions E 
of Eq. (17), that is, on the energy spectrum. The energy 
eigenvalues are the values of E that satisfy (17). The 
problem can also be discussed in terms of the dimension- 
less scaled energy lo. In this view, the solutions are the 
values of uj that satisfy (18). Using either of these two 
equations, it is easy to solve the energy spectrum numer- 
ically for a wide range of parameters. From this point 
onward, unless otherwise stated or implied, we employ 
the dimensionlcss scaled energy u. 

A standard discussion ensues from Eq. (18) when we 
plot f{uj,L) as a function of oj. In Fig. 1 the lattice 
has 16 sites, L = 16, and the added horizontal line 
depicts the value of the right-hand side of Eq. (18) for 
K, — —i. It is clear that f{uj,L) may attain any finite 
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value between two successive arguments w = — cos q for 
the given discrete set of quasimomenta q. Vertical asymp- 
totes (dashed lines) are plotted in Fig. 1 at each value of 
Lj = — COS q. For L = 16 there are eight regions bounded 
by such asymptotes, obviously because cos(g) = cos(— g). 
This docs not signal doubly degenerate energies oj, but 
is one manifestation of the fact that we have counted the 
basis vectors with q ^ and q ^ tt twice. 

For a noninteracting system with /C — )■ 0, the left-hand 
side of Eq. (18) approaches plus or minus infinity de- 
pending on the sign of the interaction parameter U, and 
the spectrum of the scaled energies is then precisely the 
numbers uj = — cosg. This is the finite-lattice analog of 
the continuum of the relative motion of two noninteract- 
ing atoms. The essential twist is that the width of the 
continuum band of energies AE = 2Mlp depends on the 
center-of-mass momentum of the atom pair P. 

Let us next restore the atom-atom interactions so that 
1/|/C| < oo. The scaled energies w that resided between 
the values of — cosg, in the band with — 1 < w < 1, will 
stay this way. These states are the finite-lattice analog 
of the dissociation continuum of the lattice dimcr. For 
brevity, we call the states with —1 < uj < 1 continuum 
states, even though there are never true continuum states 
in any finite-length lattice. The dimensionlcss energy 
values for the continuum states are labeled Wc. 

Moreover, as may be seen from Fig. 1, there is one 
energy value that moves away from the dissociation con- 
tinuum when |/C| increases and therefore 1/|/C| decreases. 
This is obviously the bound state of the lattice molecule. 
The bound state is characterized by |a;| > 1, and we 
label the corresponding dimensionlcss energy value LOb- 
The unsealed energy of the bound state E = Ml puJi, also 
varies with the center-of-mass motion, another feature 
that does not exist in free space. For attractive interac- 
tions, U < 0, the bound state is the lowest-energy state 
of the lattice dimcr and lies below the dissociation con- 
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FIG. 1. (Color online) The function /(w,L) [Eq. (18)] for 
L — 16. The horizontal line represents the right-hand side 
of Eq. (18) for = — i. The dashed vertical lines are the 
asymptotes of /(cj,L) at the values of uj = — cosg such that 
f{uj,L) = ±oo. 



tinuum. However, for repulsive interactions, U > 0, the 
bound molecular state lies above the dissociation contin- 
uum, a situation that also has no analog for molecules 
in free space. The bound state above the continuum is 
a manifestation of the trivial symmetry of the Hamilto- 
nian, H{J,U) = —H{ — J,—U), combined with the ob- 
servation that wc have already made that the transfor- 
mation J ^ — J alone has no effect on the spectrum of 
energy eigenstates. 

Having discussed the energy spectrum, we next deter- 
mine the state vectors. The generic form is specified by 
Eq. (9). The coefficients A{q) are given by Eq. (16), 
where the numerator is just a constant and will even- 
tually be absorbed into the normalization. There is an- 
other subtlety here. Namely, in quantum mechanics one 
wants to do explicit calculations using expansion coeffi- 
cients with respect to an orthonormal basis, but in the 
t J 



present case the states c\ 



c\ |0) are not orthonor- 

2- +1 2^-1 

mal for all q. The states corresponding to q and —q are 
the same, and moreover the states with q = Q and q = tt, 
which occur for integer quasimomenta, are normalized to 
two. not one. All of this requires a careful documentation 
of the range of the summation index q and, depending on 
the range used in the calculations, may lead to peculiar 
factors in the expressions governing the normalization of 
the state. The true dimension of the state space is L/2+1 
if we deal with integer quasimomenta and L/2 for half- 
integer quasimomenta. 

In our subsequent calculations wc continue the double 
counting but normalize the expansion coefficients so that 
the underlying quantum states are always normalized in 
the standard way. Now, the inner product of two states 
of the form (9) with the expansion coefficients A{q) and 
B{q) as inherited from quantum mechanics is 

(^a| Vb)=4 ^ A*(g)S((z) + 2A*(0)B(0) + 2A*(7r)B(^) 

0<g<7r 

= 2j2A*iq)B(q); (19) 

— 7r<q<7r 

if g is a half-integer quasimomentum, then the last two 
terms in the second form are missing, but the final result 
again holds true. We therefore define the inner product 
for the expansion coefficients. 



{A,B)^2Y,A*{q)B{q), 



(20) 



-7T<q<7T 



and normalize accordingly. Orthonormality with respect 
to the inner product (20) is equivalent to orthonormality 
of the underlying quantum states. Since the version of 
the Schrodinger equation (12) is Hermitian with respect 
to this inner product, the true quantum states come out 
with the proper orthonormality properties as well. 

In view of the eigenvalue equation (16) and the defini- 
tion of the dimensionlcss variables (15), the eigenstate for 
the eigenvalue to is defined by the expansion coefficients 

Ai.^q)^^^. (21) 
u! + cos q 
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Moreover, the choice of the overall numerical factor 

-1/2 



E 



{lu + cosq)'^ 



(22) 



ensures unit normalization with respect to the inner 
product (20). 

What these results mean in terms of site occupation 
numbers is seen by converting the momentum represen- 
tation state (9) back to lattice site (position) representa- 
tion using Eq. (2). The result is 



e2 



q,k\ ,k2 



(23) 

For P = this state remains unchanged in lattice trans- 
lations; for P 7^ the phase of the state changes by 
e^^ when one moves the reference point fc = back by 
one step in the lattice, meaning that fci — >■ fci -I- 1 and 
^2 — >■ A:2 -l- 1. The quantum mechanical stationary state 
of a molecule is evenly spread along the entire lattice. 
Lest this appear odd, the same applies to any aggregate 
of atoms also in free space. How the spread-out station- 
ary states relate to our intuition that a molecule is a 
localized object will be demonstrated in Sec. VB. 

At this point it is tempting to introduce the analogs of 
the center-of-mass position and the relative position for 
the two atoms, expressed in terms of the atomic positions 
ki and ^2 as 



K=\{ki+k2), fc = fei 



(24) 



Formally, we write from Eq. (23) 




(25) 



However, this expression is incorrect, formal only, since 
the summation indices K and k do not decouple in a 
simple way. For instance, K zt should be a legitimate 
lattice index, an integer. Equation (25) appears to rep- 
resent a superposition of copies of the molecule with the 
center of mass fixed at A' = 0, 



E 



a\a\\Q), (26) 



r^k 



translated along the lattice to each site A', although prob- 
lems with the summation indices persist. 

This unsuccessful attempt is as close as we have got 
to separating the center-of-mass and relative motions in 
second-quantized position representation. Nonetheless, 
we qualitatively think of the expression inside the square 
brackets in Eq. (26) as the wave function of the relative 
motion of the atom pair. 



III. LIMITING CASE OF A LONG LATTICE 

While it is easy to solve both the energy spectrum and 
the state vectors numerically for a wide range of param- 
eters, some analytical work is also possible in the limit of 
a large number of lattice sites. In the limit L — > cx) the 
quasimomenta q make an infinitely dense set over the in- 
terval (— TT, tt]. Our main technical tool is the continuum 
approximation for the quasimomenta 



E/('^) 



L 
2^ 



(27) 



valid for "sufficiently smooth" functions of quasimomen- 
tum f{q). 

Now consider the would-be continuum states in the 
limit of an infinite number of lattice sites. A simple ar- 
gument can be found by examining Fig. 1. As L — > oo, 
the curves in the domain — 1 < w < 1 get closer and 
closer together in such a way that they in effect become 
straight vertical lines. Thus, the solutions for the con- 
tinuum energies in the limit of a large number of lattice 
sites are given by the locations of the asymptotes that 
occur at 



"{Q) = ^ cos 



2t:Q 



(28) 



As before, we use the subscript c as in Uc to mark a 
consideration that is specific to the continuum states with 
— 1 < < 1, and likewise b for the bound state with 
l^bl > 1- 

The continuum energies become infinitely dense with 
L — >■ oo. An attempt to isolate any individual energy 
and eigenstate eventually becomes futile, and all observ- 
able quantities must be expressible as sums over the en- 
ergies ujc- Originally, we had the integer Q run over the 
interval (— L/2,i/2], but, as already noted, this range 
duplicates all but at most one energy eigenvalue. Given 
a "sufficiently smooth" function of the energy, g{uj), we 
therefore approximate 



L/2 

E-9(^<=) - E^ 

tJc Q=l 

rL/2 



COS 



2ttQ 



dQg 



duJr 



2ti 



This expression identifies 



27rQ 



L 



2^yr 



(29) 



(30) 



as the density of continuum states, normalized to the 
dimensionless energy of the variable lOc- 

Regardless of the ultimate limit of an infinitely long 
lattice, L — > oo, we always think of the lattice as finite. 
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All inner products, normalizations, and so on, are with 
respect to discrete sums. When appropriate, these sums 
are just approximated as in Eqs. (27) and (29). 



A. Bound state 

First we will analyze the energy and the normalization 
of the bound state. In the limit of a large number of 
lattice sites, we approximate 



1 



uJb + cos q 



L 1 
2^L 



1 



ujb + cos q 



(31) 



a valid process since for any /C 7^ the bound state has 
\u!fj\ > 1 and the function of q to be summed does not 
become singular in the limit L — > 00. We therefore find 
the bound-state energy by substituting (31) in (18) and 
solving for ujb, 



Wb = sgn(/C)v/l+/C2 



(32) 



In the limit of strong atom-atom interactions, this re- 
duces to w;, ~ /C and so in dimensional units the energy 
of the bound state is E), = hU, but while |/C| ~ 1 the 
energy of the bound state depends on the center-of-mass 
motion. These results for the bound-state energy agree 
with those found, for example, in [8] and [14], although 
these authors do not mention the limit L — > 00. 

The state vector for the bound state is of the form 
of (21) and (22) with cj = Wf, = ±Vl + IC^- To find 
the normalization coefficient in the limit L — > cx), we 
calculate 



[CbiuJb)] 



2L\LJb\ 



(33) 



The bound state is specified by (21), (33), and (32) as 

|/C|3/2 1 



Ab{ujb,q) 



^j2L\uJb\ + cos q 



(34) 



B. Continuum states 

We next consider the state vectors for the continuum 
states. While the bound state has a straightforward limit 
as L — >■ c», the continuum states with — 1 < Wc < 1 pose 
a problem. For any finite number of lattice sites the 
stationary state is, of course, still of the form 



Ac{u!c,q) 



Wc -I- COS q 



(35) 



However, as the continuum states get denser with L — >■ 
00, this expression becomes singular as a function of q 



in the neighborhood of u>c + cosg ~ 0, and a rule for 
handling the singularity has to be defined. 

Our assignment is basically the following: In principle 
we only do discrete sums, and knowing the exact dis- 
crete continuum eigenvalues ujc would always allow us 
to carry out the sums without problems. Nevertheless, 
in practice we would like to approximate the sums as 
integrals, as in (27). We arc therefore looking for a func- 
tion of (7, Ac{u)c,q), which would give the same value 
from a continuum-limit integral (27) as does a discrete 
sum involving the amplitudes (35). We will, in fact, find 
such a function Ac{u!ct<i) [Eq. (50)]. Moreover, it comes 
with the appealing property that, unlike in applications 
of Eq. (35), one does not have to know the energy uJc 
precisely. 

Following Ref. [19], we try the following ansatz: 



1 



= P- 

uJc + COS q tjJc + cos q 



+ Tl{ujc)5{uJc + cos q) , (36) 



where P stands for principal- value integral and Il{u>c) is a 
function to be determined. To find 11 (wc), we attempt to 
solve the continuum version of the Schrodinger equation 
itself. By rearranging Eq. (16) and using the dimension- 
less variables (15), it is 



{uJc + cos q)Ac{q, UJc) ^ — J dq' Ac{uJc,q') ■ 



(37) 



Substituting Eq. (35) into (37) and taking into account 
that for UJc G (—1,1) we have 



[uic + cos q) P = 1 , 

UJc + cos q 

{uJc + cos q) 5{uJc + cos q) ~ 0, 

' 0, 



dq ■ 

UJc + cos q 

we immediately find that 



n(c^e) 



K 



(38) 
(39) 
(40) 

(41) 



leads to the eigenvalue uJc in Eq. (37). 

Next wc digress on the normalization of the contin- 
uum wave functions. In the original discrete case two 
state vectors corresponding to two different energies are 
orthonormal in the sense of a Kronecker 5, 



2 ^ Aciujc, q)Ac{uj'c, q) = 5^^^^^'^ 



(42) 



Taking an "arbitrary" function F{ujc), wc therefore have 
^ ^ Aci^c. q)Ac{uj',, q)FU) = \ Fiujc) . (43) 



As one more remnant of the state counting problems, 
there are only half as many energies uJc as there are co- 
efficients Ac{uJc, q)- The sum over uj'^ is to be understood 
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accordingly. On the other hand, we have 

EL dujr 
2n Jl-uj^ 



(44) 



Using the continuum approximation for the sum over q 
in Eq. (43) once more, we have 



1 dc^^FK) r 

-1 



dq Ac{uJc,q)Aciuj'^,q) 



The correct continuum approximation normahzation of 
the coefficients Ac{qTLOc) therefore reads 



dqAc{LOe,,q)Ac{Lo'^,q) 



(46) 



To implement this normalization, we start from 
Eqs. (35) and (36), substitute x ~ —cosq, and find 



dqAc{uJc,q)Ac{ui'^,q) = 2Cc(wc)Cc(Wc) 



dx 



-Il{uic)S{u;c-x) 



■ + Uiu;',)SU~x) 



(47) 



By virtue of the form of the function n(a;c), Eq. (41), the integrals of the products involving a principal value and 
a 6 function cancel, and the integral of the product of two 5 functions is simple. The product of two principal- value 
integrals can be handled with the identity [19] 



1 



UJr ~ X LO' — X UJr — LO' 



1 



P- 



1 



+ 7r^(5(x — ujc)5{x — u}'^) . 



The first two terms both give zero by virtue of Eq. (40), so that we are left with 



dqA,{u:,,q)A,U,q)=2CciuJc)^ + n(^c)' 



(48) 



(49) 



Comparison with Eq. (46) gives the normalization coefficient Cd^c), and hence the complete state vector for the 
continuum state 



Ac{ujc,q) = 



> 1 

OJc + COS q 



'-5{lOc + cosg) 



(50) 



The precise meaning of Eq. (50) should be borne in 
mind very carefully. This is an approximation to the 
continuum wave function to be used when one replaces 
the sum over discrete quasimomenta with the continuum- 
limit integral as in Eq. (27) and is purportedly valid in 
the limit of a large number of lattice sites, L — > oo. At 
this point one, in principle, needs to know the continuum 
eigenvalue Uc- However, if the final result of the cal- 
culation varies smoothly with Wc, one may equally well 
think of Wc as a continuous energy with the density of 
states (30). 

IV. DIMER IN POSITION REPRESENTATION 

For comparison and additional results, we next study 
the lattice dimer in position representation, somewhat 
parallel to Ref. [14]. The starting point is the time- 
independent Schrodinger equation (12) in the case when 
both and q are legal quasimomenta. Some changes 
in the formulation would result for half-integer quasimo- 



menta, but we do not discuss them. 

As already noted in connection with the somewhat ill- 
defined equations (25) and (26), we use 

«fe = ^E^"''^(9) (51) 

to represent the relative motion of the two atoms in posi- 
tion representation, k being the distance in lattice units 
between the atoms. From another viewpoint, the step 
from A(q) to ak in Eq. (51) is a discrete Fourier transfor- 
mation, a perfectly well-defined mathematical operation. 
The transformation preserves the inner product (20), so 
that the position-representation states are normalized ex- 
actly as in momentum representation. 

Transforming Eq. (12), we directly have an equation 
for the amplitudes ak, 

ojak + ^ak+i + \ak^i = 4,o A^ao • (52) 

Here we choose the values of the relative coordinate of the 
two atoms from the interval k € (— L/2,L/2], but more 
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generally, ak must be regarded as a periodic function of 
k with the period L. 

Now, Eq. (52) is a second-order finite-difference equa- 
tion that may be solved like the corresponding second- 
order differential equation. The homogeneous equation 



ujak + hak+i + hak-i = 



(53) 



admits solutions of the form = x'^ , where we see from 
direct substitution that the constant x may take on the 
values 



(54) 



The product of the two possible values of x is always 
unity. 



A. Bound state 

First consider the case \lu\ > 1. In this case, lu ^ ui, 
and X is a real number. Suppose > 1, then the two 
values of x smaller and larger than unity in absolute value 
are 



x< = -Wfe + - 1, Xy = -Ub - yuj^ - I ■ (55) 

By the boson symmetry, we may always require that 
= Q!-fc. Recalling that a;> — x^^, the only possible 
solution to Eq. (52) is then of the form 



ak = 



Ax'^ + Bxi^ , fc > , 
Ax^ -t- Bx%^ , fc < . 



Substitution into Eq. (52) with fc = gives 



{A-B)JujI-1 = {A + B)K.. 



(56) 



(57) 



However, there is another equation to reckon with. 
Namely, to satisfy the periodic boundary conditions, the 
solution must be periodic, so that at fc = L/2 the ex- 
pression (56) must switch between the fc > and fc < 
forms while remaining a solution to (53). This leads to a 
second condition, 



(Ax</^ + Bxy''^^ + \ (AxI 



L/2+1 



Bx. 



-L/2+1 



~Ax 



L/2-1 



Bx 



L/2-1 
> 



B = Axb . 



(58) 



(59) 



In the limit of an infinitely long lattice, L — > oo, we have 
B/A = 0, so that Eq. (57) gives 



Wb = sgn(/C)Vl+'*C2 



(60) 



Here we have worked out the case /C < as well, which 
leads to Wf, < —1. The result, of course, is as before 
[Eq. (32)]. The state vector for the bound state is [14] 



akiujb) = Ax 



\k\ 



JC 

2a;f, 



(/C - uJb) 



\k\ 



(61) 



where we have explicitly normalized, ^J^k lo^fcP = 1- 

If the lattice is not infinitely long, the bound state also 
has a component that grows exponentially away from the 
center site fc = as |a;>|l*''l. Moreover, the bound-state 
energy is shifted from the infinite-lattice value. In fact, 
Eqs. (57) and (59) have a solution for A and B if and 
only if 



1-IC 



1 + /C) 



(62) 



If this equation for Wf, has an analytic solution for a gen- 
eral (even) L, we have not been able to find it. However, 
the equation poses no particular problems numerically, 
and we have verified that the solution agrees with what 
is found from Eq. (18). The coefficients A and B can be 
determined explicitly in terms of oj;, so that the state (56) 
is properly normalized, but the best we have been able to 
do analytically is cumbersome. Generally speaking, if a 
finite lattice is the issue, we believe that there is little to 
be gained from an attempt to press analytical as opposed 
to numerical calculations. 



B. Continuum states 

Second, consider the case |a;| < 1, so that uj = uJc- 
Then x< and x> are complex numbers of unit modulus, 
of the form e^*'' for some real q; 



(63) 



The solutions to Eq. (52) with the requisite boson sym- 
metry are 

, ■. _ ( Acoskq + B sinkq, fc > 0; 
'^'^y^'^) - \^ A cos kq- B sin kq, fc<0. ' 

From Eqs. (63), the corresponding energy is of the form 
LOc ~ ~ COS q . (65) 

Equation (52) with fc = then gives 

Bsinq^lCA. (66) 

Once more, the ansatz (64) must give a solution to 
Eq. (53) also at fc = L/2, where the solution switches 
between the two forms, so that we have 

UJc {A cos [(iL) q] +Bsm [{^L) q] } 

-|-{Acos [(ii - 1) q] +Bsm [(iL - l) } = , (67) 
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h(q, L) 

3 r 



FIG. 2. (Color online) The function h{q,L) [Eq. (69)] for 
L = 16. The vertical dashed lines represent the asymptotes 
of h{q, L); the horizontal line stands for IC — ^. 



or 



B = Atan^. (68) 

The quantization condition for q now becomes 

/i(g, L) = sin q tan ^ = K. . (69) 

We plot the function h{q,L) for L = 16 in Fig. 2; the 
horizontal line depicts the value of the right-hand side of 
Eq. (69) for /C = |. The plot shows that the function 
h{q, L) attains every real value > precisely L times 
when q varies in the interval (— tt, tt), once in each interval 
of the form {ttQ/L, ttQ/L + tt/L) for an integer Q. 

Let us again for definiteness take /C > and consider 
nonnegative roots q G [0,7r) to Eq. (69) (if g is a root, 
then so is —q). Such roots for /C = would be 



PQ 



2ttQ 
L ' 



Q = 0,...,i/2-l 



(70) 



We characterize the roots qg for /C > with the same 
"quantum number" Q in such a way that qg develops 
continuously from pg when /C is continuously increased 
from zero. In the limit L — > oo the change in sin q is neg- 
ligible for the change of K, from zero to its final value, so 
that the actual energy value satisfies uJciQ) = — cosqg ~ 
— cos pq, and therefore sin qq ~ -^/l — cos^ pq. Noting 
this, one immediately sees that to the leading nontrivial 
order in 1/L, the Qth root of Eq. (69) is 



Qq -Pq + Y arctan —= 



COs2 Pq 



PQ 



(71) 



The branch of the explicit arctan function is chosen so 
that the value lies between — ^tt and ^tt, and the addition 
of PQ corresponds to the choice of the proper branch of 
the arctan function to solve Eq. (69). 

After this organizational work, the roots of Eq. (69) 
and the corresponding energies uJciQ) — — cosqg are 



easy to find numerically. The results agree with those 
obtained by solving Eq. (18) numerically. Again, if the 
finite number of lattice sites is the issue, we recommend 
direct numerical computations. 

The state vector in momentum representation is the 
inverse of the discrete Fourier transformation (51), or 



1 

7l 



L/2-1 

E ' 

fc=-L/2 



''P^a.iiOc), (72) 



where the momenta p, in fact, run over the values called 
PQ in Eq. (70), albeit with Q = -L/2 +1,..., L/2. In- 
serting Eq. (64), making use of Eqs. (68) and (69), and 
noting that sin(pL/2) = 0, we have 

A,{uj,,p)^^ , (73) 

VL cosp — cosg 

where D is a so-far undetermined normalization constant. 
Given that cosg = —coc, this is in agreement with the 
previous expression for the state vector (35). 

The problem in an attempt to take the limit X — > oo 
is the terms with p q, which make the expression (73) 
singular. However, by virtue of Eq. (71), we are now in 
a position to take the continuum limit differently than 
we did in the momentum representation. Following an 
old [20, 21] but maybe not so generally known idea, let 
us consider Eq. (73) in the vicinity of the singularity that 
occurs around a given q value labeled qg^ , with p > and 
q > for definiteness. We then have 

A,{wc,pg)^ y=-. (74) 

VL sinpQ„{pg - pg„ - AqJ 

In the limit L — !■ oo, an arbitrary fixed small symmetric 
neighborhood Ap of the momentum pg^ contains a very 
large number of momentum modes pg. Taking a smooth 
function of momentum G{p)^ we estimate 

^ Ac{uJc,Pq)G(pq) 

PqEAp 



DfC 27r 
VX isinpQo 

DK 27r 
y/Z L sin pQg 



G(pQo)E 

Q= — co 



Q-Qo-LAQj{27r] 

^AQ.,/(27r) 
.Q2_(iAQj(27r))2 



(75) 

from duplication of the original sum with the change of 
the summation index Q — > —Q. However, the sum [22] 



E 

k—'-OO 



-TT cot[7ra;] 



(76) 



and the definition of Ag„ from Eq. (71) give 
Ac{ujc,pg)G{pg) 



DK. L 



G{pQo) 



1 TT^l - COsVqo 



(77) 
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On the other hand, replacing the sum with momentum 
points spaced at the intervals 27r/L by an integral, we find 



27r 

T 



PIC 



1 TT^A 



JC 



G{PQo)- 



(78) 



Clearly, in the small range Ap, the function Ac{uJc,p) 
behaves as 



Ac{u)c,p) 



DK 



1 TT^ 1 - C0S2 PQ^ 



7l 



JC 



■ S{cOSp + Wc) 



(79) 



where we have noted in the last step that Wc — — cos pq^ . 
Other than in the small and symmetric neighborhood 
Ap, Eq. (73) still applies. In the limit L — > oo, the 
momentum representation state vector therefore behaves 
under integrals over p as 



Ac{ujc,p) 



PIC 

7l 



1 



ijJc + cos p 



■ ^(cJc+COSp) 



(80) 

The result has the same functional form as before 
[Eq. (50)] so all that remains is to verify the normaliza- 
tion constant. To this end we first note that the condition 
for normalization of the state (64) 




L{K? + 1 - ^2) _ 2KuJc 



(81) 



(82) 



by using Eqs. (68) and (69). In the limit L ^ oo, Eq. (50) 
immediately follows. 



V. EXAMPLES OF DIMER DETECTION 

A. Momentum distribution 

In our first foray into the detection of the dimers we 
assume that after a preparation of a possibly large num- 
ber of bound-state dimers, the lattice is removed and 
the atomic cloud expands ballistically. The positions 
of the atoms are then detected after some free-flight 
time. Ideally, this procedure converts the momentum 
distribution of the atoms into a position distribution so 
that a measurement of the position distribution amounts 
to a measurement of the momentum distribution. Our 
thought experiments closely mimic actual laboratory ex- 
periments [8]. 



There is a complication arising from the periodicity of 
the lattice that was also discussed in Ref. [8]. Namely, 
if the lattice is switched off instantaneously, the momen- 
tum distribution consists of periodic repetitions of the 
first Brillouin zone modulated by the momentum distri- 
bution of an atom in one individual lattice site, that is, 
the momentum distribution of the one-atom states asso- 
ciated with the annihilation operators ak- In the experi- 
ments [8] the lattice was turned off on a time scale such 
that the structure of the physics on a length scale below 
one lattice spacing was presumably removed adiabati- 
cally, while the site-to-site physics did not have time to 
adjust. The result is a momentum distribution confined 
to the first Brillouin zone. We analyze such distributions 
as well. 

Now take an eigenstate of the center-of-mass mo- 
tion (9) with the total momentum P. A straightforward 
exercise then gives the momentum distribution. 



Mip) = (^1 clcp \^) = 4 \Ab {ujb,p- hP) I 



(83) 



with the A function given in Eq. (34) for the bound state. 

First consider a stationary center of mass, P = 0; then 
we have 



M{p) cx 



1 



(cjfc + cospY 



(84) 



For a repulsively bound state, > 1 holds true, so that 
the momentum distribution is a maximum at the edges of 
the first Brillouin zone, when cosp = —1, ov p = ±7r. For 
the bound state with attractive interactions, the maxi- 
mum is at the center, p = 0, of the first Brillouin zone. 
These features were seen experimentally [8]. 

At first it might seem that the variation of Eq. (83) 
with the center-of-mass momentum P would simply be 
to slide around the momentum distribution of the atoms 
cyclically in the first Brillouin zone by \P. This is not 
the case, since the energy of the bound state Wb also de- 
pends implicitly on P. For convenience, we define a mo- 
mentum distribution /(p, P) = LAAiji) / At: normalized so 
that J^^f{p, P) dp = 1 and have the explicit expression 



fip.P) 



|/C(P)|= 



2TT^T+iap)^ 
1 



{cos(p - iP) + sgn[/C(P)]yTTX(Pf }■ 



, (85) 



with 



/C(P) = 



/Co 



cos ^P 



U_ 
2J 



(86) 



We plot an example with /Co = 8 in Fig. 3. This is a con- 
tour plot with p as the horizontal axis, P as the vertical 
axis, and brighter shades standing for larger values. The 
ip sliding of the distribution of the atomic momenta p 
with the center-of-mass momentum P is visible, but there 
is also a modulation so that, if anywhere, the momentum 
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-1.0 -0.5 0.0 0.5 1.0 

p/jT 



FIG. 3. (Color online) Contour plot of the momentum dis- 
tribution in the bound state of the lattice dimer, /(p, P) of 
Eq. (85), as a function of the momenta of the individual 
atoms p and the center-of-mass momentum P, for U/J = 16. 
Brighter shades represent larger values. 



distribution of the atoms is always peaked near the edges 
of the Brillouin zone. The peaks in the momentum dis- 
tribution are narrowest for the center-of-mass momenta 
P = — 27r, 0, and P = 2tt; at P ~ ±tt, the momentum 
distribution is completely flat, /(p, ±7r) = l/27r. 

A discussion of the momentum distribution as a func- 
tion of the center-of-mass momentum based on numerical 
solutions of the lattice as a many-body system was also 
offered in Ref. [8]. It brings up similar qualitative ele- 
ments as our discussion, but it seems to us that the trans- 
lation of the momentum distribution with the center-of- 
mass momentum is quoted in Ref. [8] as P, whereas we 
obtain iP. 



B. Size of the bound state 

Our next example on the detection of the dimer is 
a thought experiment in which there are precisely two 
atoms in the lattice, and the number of the atoms at 
each site is measured. This experiment is carried out 
over and over again, and the detection statistics is com- 
piled. In our example we assume an absorbing detector 
that removes an atom from further consideration once it 
has been observed. Modeling after the well-known pho- 
ton detection theory [23, 24], the joint probability for 
finding an atom at sites ki and k2 in the energy eigen- 
state of the lattice dimer (23) with a fixed center-of-mass 
momentum P is 

J(fci,fc2) = N {u}\ al^al^a^^ak, \uj) 



= — (87) 

where is a normalization constant. As befits trans- 
lational invariance, the probability only depends on the 
distance between the lattice sites. The distribution of the 
distance is governed by the internal state of the molecule 
in position representation, ak of Sec. IV. 

We will not analyze the implications for a finite size 
lattice, but go directly to the limit of infinite lattice with 
L — >■ oo. By virtue of Eq. (61) the variance of the distance 
between the detected atoms in the bound state is 

Any eigcnstate (25) for a fixed center-of-mass momen- 
tum P of the two bosons in the lattice is translation in- 
variant. Taken individually, the detected atoms must be 
evenly distributed along the lattice, and there is no sign 
of the internal state of the dimer. The internal state 
is manifested in the correlations between the detected 
positions of the atoms; in the ground state the atoms 
are observed in pairs with the root mean square distance 
l/(-\/2|/C|) between them. Even if the interactions be- 
tween the atoms are repulsive, the stronger the interac- 
tions, the more tightly the atoms are paired. 



C. Dissociation rate of the bound state 

In order to find the bound-state energy experimentally, 
some form of spectroscopy has to be carried out. Here we 
assume that the intensity of the lattice light is modulated 
periodically. The dissociation rate of a repulsively bound 
pair (U > 0) as a function of the modulation frequency 
was studied experimentally in this way in Ref. [8]. 

Given that the tunneling amplitude is much more sen- 
sitive to the depth of the optical lattice than the atom- 
atom interactions, in our model only the hopping matrix 
element J becomes time dependent: 

J J{1 + Acosvt) , (89) 

where A J and v are the amplitude and frequency of the 
modulation, respectively. The modulation adds a "per- 
turbation" to the Hamiltonian. In lattice momentum 
representation, we write it as 

H' 

— = Acosi/t ^Wqcjcg . (90) 

Q 

The matrix element of this perturbation between any 
two states i and j of the form (9) turns out to be nonzero 
only if the states have the same center-of-mass momen- 
tum P; then we have 

H' 

(V-il IV^j) = -AQpMij cosvt . (91) 



12 



It proves expedient to define dimensionless matrix ele- 
ments, scaled to the frequency flp, as 



Mrj =2j2cosqA*{q)A,iq) 



dq cosqA*{q)Aj{q) 



(92) 



Some results in what follows have maybe unexpected 
complications in the notation, which are needed to cover 
the possibility that Hp < 0. 

The relevant L — )- oo energy eigenstates are specified 
by Eqs. (34) and (50) . The matrix elements are 



M, 



bb 



1 

UJb 



M, 



be 



(93) 
(94) 



Mr; 



L 



5{uJc - UJc') 




(l-c.,)(l-^0 



(/C2 + l-W,)(/C2 + l-W^) 



(95) 



The matrix elements Mcc' may look severely singular, but 
they are not. Namely, recognizing 27r^l — uj'^/L as the 
inverse of the density of the states Wc, the finite-lattice 
version of the delta function part is 



(96) 



This is a perfectly well behaved diagonal matrix element 
for the continuum states. 

The form for the amplitudes A in Eq. (50) was sup- 
posed to be good inside integrals only if the rest of the 
integrand is well behaved, which is clearly not the case 
in Eq. (92) when i and j both stand for continuum 
states. The matrix elements (95), however, transcend 
their derivation, with the understanding (96) that they 
agree with the matrix elements computed numerically in 
a large finite lattice. 

Writing the state |w) as \b) for the bound state and as 
|c) for the quasicontinuum state, the Hamiltonian written 
in the eigenbasis of is 



H 



: UJb 



-A cos vt 



M, 



bb 



J2^c\c){c\ 



Mcblc 



(97) 



We employ perturbation theory in the dimensionless pa- 
rameter characterizing the modulation depth. A, to study 
the dissociation rate of the bound state to the quasicon- 
tinuum states. Here the problem of time-dependent per- 
turbation theory is unusual in that there are diagonal 



transition matrix elements. We attempt to get past this 
obstacle with the assumption that the system starts in 
the bound state. The matrix elements Mbc and Mcb must 
then be involved. The corresponding terms in the Hamil- 
tonian are already in themselves first order in the small 
parameter A. We therefore apply the Hamiltonian in the 
form 



H 



= ojb\b){b\+J2'^c\c){c\ 

C 

-A cos iytJ2iMbc\b){c\ + Mcb\c){b\) (98) 



in the hope that we get a correct description to lead- 
ing order in perturbation theory in the parameter A. 
This should be warranted if the coupling between the 
continuum states mediated by the matrix elements Med 
does not cause the exact solution of the time-dependent 
Schrodinger equation to become a nonanalytic function 
of the parameter A at A = 0. We assume so without 
further ado. 

At this point the standard Golden Rule transition rate 
applies and gives the dissociation rate of the bound state 
in proper dimensional units as 



r = ^^^MM,(A). 



(99) 



Here the energy-conserving continuum state with the la- 
bel Wc = A depends on whether the bound state lies 
above (— ) or below (-f ) the continuum. 



p rip 



(100) 



This is an analog of a parameter called detuning in laser 
spectroscopy, and would be controlled in practice by 
varying the modulation frequency v. Using the density 
of states given in Eq. (30), we have 



r ^ 



\Qp\\JC\^VT^l 
2K|(c.2-A2) 



(101) 



The shape of the dissociation line, variation of the dis- 
sociation rate with the modulation frequency, changes 
from (1 — A^)^^/^ in the limit of weak interactions 
|/C| < 1 to (1 - A2)i/2 in the limit |/C| > 1. In di- 
mensional units the width of the spectrum is 2|f2p|. The 
total strength of the transition is characterized by 



dA 



r(A) 

" \ 2(1 + IC^){\IC\ + VTTk^) 



(102) 



which tends to 7r|/Cp/2 as |/C| < 1 and to 7r/4 as |/C| > 1. 

The experimental results [8] are mostly in the limit 
/C ^ 1. While there is no indication in Ref. [8] whether 
perturbation theory in the modulation amplitude should 
apply or not, our results qualitatively explain the pub- 
lished line shape. 



13 



VI. CONCLUDING REMARKS 

We have demonstrated how various molecular physics, 
condensed-matter physics, and many-body physics as- 
pects come together in a description of a dimer of two 
bosonic atoms in an optical lattice. Aside from a sur- 
prising number of mathematical complications that we 
have sorted out both explicitly and behind the scenes, 
the main technical issue here is the old [20, 21] quasi- 
continuum problem: what to do with a system that has 
a dense set of energy eigenstates? We have described 
the stationary states of the dimer and found analytical 
expressions for the stationary states in the limit of an 
infinitely long lattice. Once the groundwork is laid, ap- 
plications are straightforward. As an example, we have 
briefly discussed three ways to detect a bound dimer. 

We have analyzed the one-dimensional Bosc-Hubbard 
model as a specific example. Nevertheless, we can think 
up, and several authors [9, 10, 15, 25] have thought up, a 
number of similar problems. Ultimately, what we hope to 
have achieved here is a template, a methodology, that ap- 

I 



plies to all sorts of dimer problems in lattices and will pos- 
sibly contribute to future studies of aggregates of atoms 
in a lattice with more than two atoms as well. 
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Appendix A: Operator Algebra Example 

The second-quantized notation automatically takes 
care of the symmetries with respect to particle exchange, 
but we may run into possibly long products of creation 
and annihilation operators. In some subfields of physics 
the methods for dealing with them are standard fare, but 
for completeness we demonstrate a technique employed 
frequently in this article by deriving Eq. (8). First, given 
Eq. (7), we have from the left-hand side of Eq. (8) 



E ^1 



91-1-92,93+94] ''iji '^92*^93 '^94 



IV') 



91,92,93,94 



E 

91 ,92, 93, 94, Pi ,P2 



91+92,93+94] '^9l''92''93'^94Cpi'-p2 



\cl 10) 



(Al) 



The general idea is to put the products of the operators into normal order, annihilation operators to the right of 
creation operators, by using the commutators (or anticommutators) of the operators. Here CpcJ = cjcp -I- S[p^q] holds 
for arbitrary p and q. We therefore have the following chain of manipulations: 



^\l^\-2'^1i^<li^'pi^'p2 |0) ~ ^\-L^\2^^<13^'pi^1i^'p2 '^[94,Pl]''93''p2) 1^) ~ "^ijl '^L ('^[94 ,P2] "^93 "^jji ^ '^[94 ,Pl ] '^[93 ,P2] ) 1^) 



94,P2l"l93,PlJ 



■'5[94,pi]'5[93,P2]) |o) ; 



(A2) 



once an annihilation operator has been moved to operate on the vacuum ]G), the corresponding term vanishes. Sorting 
out the Kronccker J's and using the symmetry A{pi,p2) ~ A(p2,Pi), (Al) and (A2) immediately combine to give the 
right-hand side of Eq. (8) . 
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